The dynamical foundation of fractal stream chemistry: The origin of extremely long 

retention times 
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We present a physical model to explain the behavior of long-term, time series measurements of 
chloride, a natural passive tracer, in rainfall and runoff in catchments [Kirchner et al, Nature, 
403(524), 2000]. A spectral analysis of the data shows the chloride concentrations in rainfall to 
have a white noise spectrum, while in streamflow, the spectrum exhibits a fractal 1// scaling. The 
empirically derived distribution of tracer travel times h(t) follows a power-law, indicating low- 
level contaminant delivery to streams for a very long time. Our transport model is based on 
a continuous time random walk (CTRW) with an event time distribution governed by ij}(t) ~ 
A^- 1 - 13 . The CTRW using this power-law tp(t) (with < (3 < 1) is interchangeable with the time- 
fractional advection-dispersion equation (FADE) and has accounted for the universal phenomenon of 
anomalous transport in a broad range of disordered and complex systems. In the current application, 
the events can be realized as transit times on portions of the catchment network. The travel time 
distribution is the first passage time distribution F(t; I) at a distance I from a pulse input (at t = 0) 
at the origin. We show that the empirical h(t) is the catchment areal composite of F(t; I) and that 
the fractal 1/f spectral response found in many catchments is an example of the larger class of 
transport phenomena cited above. The physical basis of tp(t), which determines F(t;l), is the origin 
of the extremely long chemical retention times in catchments. 



INTRODUCTION 

The distribution of travel times required for chemi- 
cal substances to travel through a catchment is an im- 
portant hydraulic quantity, determining the retention of 
pollutants until they are eventually released to a stream 
or lake. The same distribution controls the transport of 
chemicals in subsurface hydrological systems. The na- 
ture of this distribution determines the expected ecolog- 
ical impact of contaminants. 

In a long-term catchment study conducted at Plyn- 
limon, Wales, Kirchner et al. [2000] (hereafter denoted 
KFN) compare the (input) time series of the concentra- 
tion cji(t) of chloride tracer in the rainfall to the (out- 
put) time series of the concentration c s (t) into the Hafren 
stream. They relate the concentrations through the con- 
volution integral 



spectral power of the water fluxes of rainfall and stream- 
flow coincide while the spectra of c s (t) show strong atten- 
uation on all wavelengths less than 5-10 years. Further, 
the cn(t) spectra scale as white noise in sharp contrast to 
the fractal 1/f scaling of the c s (t) spectra. Using these 
results in the Fourier transform of the convolution in (Q) , 
with / denoting frequency, and C s , H and Cr represent- 
ing the respective Fourier transforms, 

C s (f) = H(f)C R (f) (2) 

and the relation of the power spectra is 

\C s (f)\ 2 = \H(f)\ 2 \C R (f)\ 2 (3) 

Since |Cr(/)| 2 is nearly a white noise spectrum (a con- 
stant) one has \H(f)\ 2 oc |C S (/)| 2 . From the measured 
power-law scaling of \C s (f)\ 2 oc /~ ' 97 , KFN conclude 
that 



,(<)= f™h(t')c R (t-t')dt' 



(1) 



h{t) ~ r m 

)), where m 



(4) 



where the effective travel time distribution h(t) governs 
the lag time between injection of the tracer through rain- 
fall and outflow to the stream. 

KFN observe that the streamflow (volume of water) 
responds rapidly to storm rainfall inputs while the c s (t) 
response is highly damped (i.e., low-level contaminant de- 
livery to streams for a very long time) . The timescales of 
catchment hydrologic and chemical response are resolved 
using spectral methods || (in which the input /output re- 
sponse can be compared at each (time) wavelength) . The 



(from H(f) ~ f~ {1 ~ m> ), where m « 0.5. 

In order to ensure that h(t) is integrable and possesses 
a finite average travel time, KFN arbitrarily chose the 
gamma distribution 



hit) 



-t/ e 



(5) 



^r( 7 ) 

with 7 = 0.48 and g = 1.9 years, with a characteristic 
travel time Th = L h(t)tdt — jg. Due to the rela- 
tively large value of g, essentially the entire data corre- 
spond to the power-law h(t) ~ i 7-1 with the power spec- 
trum |C S (/)| 2 cx \H(f)\ 2 - /- 2 t. KFN refer to similar 
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scaling between / _0 ' 7 and f~ 12 found in Scandinavian 
and North American field sites, which indicates a certain 
ubiquity to fractal, scale-free forms (which is the issue we 
address in this letter). 



TRANSPORT MODELS 

The first step to understanding transport in the catch- 
ment is the clarification of the meaning of (|l|). The h(t) 
is the effective response to a pulse of rain falling on the 
entire area of the catchment. Every point of this area is 
a source of chloride and the stream is a line sink for the 
chloride. We simplify the area to be a rectangle of width 
A about this stream sink (with periodic boundary condi- 
tions on the sides perpendicular to the stream |L(|). The 
sink or absorbing boundary ensures the correct "count- 
ing rate" of chloride at the arriving point, i.e., the first 
passage time distribution F(t; I) is the appropriate travel 
time distribution from a pulse source at point I. In terms 
of this intrinsic distribution, 



c s (t; l a 



E 



F(t';l-l s )c R (t-t';l)dt' (6) 



where c s (t; l s ) is the chloride concentration at the stream 
position l s at time t, Q is the size of the catchment and 
cr(£; I) the rain- input at a position I in fl. We can con- 
sider the sampling positions {/ s } of c s to be a small region 
compared to fl and hence c s (t;l s ) w c s (t). To do the l- 
sum we assume that the rainfall is distributed uniformly 
in fl, Cfj(t; I) « Cr(£). Hence, we recover (Q) now with 



(7) 



The distribution F(t; I) must be integrable in time and 
thus so hit) for a finite tt. The basis of comparison for 
various transport approaches is the computation in (]?]). 



Advection-dispersion equation (ADE) 

The textbook approach to transport of passive tracers 
in both surface and subsurface systems usually focuses 
on the advection-dispersion equation (ADE) 



dW/dt + vdW/dx = Dd 2 W/dx 2 



(8) 



e.g., for Id with constant v the average fluid velocity 
and D the dispersion coefficient. This equation governs 
the temporal evolution of the probability density function 
(pdf) W(x,t) of finding a tracer particle, which under- 
goes dispersive motion, at a certain position x at time 
t after release at t = 0. The pdf W is the normalized 
concentration profile. 

The travel time distribution belonging to (^) can be 
obtained by the Laplace Transform (LT) technique. This 



yields the LT P(u;x) = T{t;x)e~ ut dt of the first 
passage time from x to the absorbing boundary at x — 0, 
T{u; x) = exp (x [v - y/ADu + v 2 \ /(2D)). LT inversion 
delivers the travel time distribution 

which when inserted for F(t; x) in (0) yields 

h(t)oc [ !F(t; x)dx = — [erf (z) + erf [ — - — z ) 
Jo 2 J 

, exp(-2: 2 )-ex P (-(^~z) 2 ) 



where the Peclet number P e = Xv/D, z = \\Jv 2 t/ D 
and erf(z) is the error function (the same result (|10|) was 
also derived in Kirchner et ai, 2001). For t < 0.1D/v 2 , 
h(t) ~ t~i which is similar to (jij) for m ps 0.5. [This 
upper limit on the time range holds for P e > 1 ; for P e < I 
the limit is smaller.] The range is estimated by assuming 
D to be the macrodispersion, D = vol, where a is the 
dispersivity. Assuming v « I0 2 m/yr and on the km scale 
a - (1 - 10 2 ) m, t < 0.1a/ v - (10~ 3 - lO" 1 ) yr. One 
needs v w 10 m/yr and a ~ 10 2 m in (10), for h(t) ~ 



to partially overlap the time range of observation of this 
dependence [Q. Moreover, this behavior is only the same 
as ([|) in the special case where m (or 7) = |. The 
catchment data (spectral power ~ / _27 ) quoted above 
cover a range of 0.7 < 27 < 1.2. 

For large t, h(t) in ( |l0| ) exhibits an exponential de- 
crease ~ t~2 exp (— v 2 {AD)~ 1 t) which is faster than (0). 
The decrease in hit) at large t it >> X/v) is due to the 
finite size of the catchment and not an arbitrary limiting 
time. The relation between the operative time range of 
(0) and A will be determined in the following. 



Continuous-time random walk and 
time- fractional- ADE 

As shown above the standard treatment of the trans- 
port using the ADE yields a marginal accounting of 
the data. A transport model must determine a F(t; I) 
that when inserted in (^) gives rise to (Q) for over three 
decades of time [0 . As we will show this is best expedited 
with a non-Gaussian form of Fit; I). A particularly dis- 
persive form of Fit; I) is associated with anomalous trans- 
port which has been observed in a wide range of disor- 
dered and complex systems: electron hopping/multiple- 
trapping in amorphous semiconductors |l(J, particle mi- 
gration in fracture networks M , contaminant transport in 
heterogeneous porous media |l|, anomalous diffusion^, 
and Hamiltonian chaos Q. There are key common fea- 
tures in these anomalous transport phenomena, e.g., non- 
Gaussian propagation of an initial pulse of particles with 
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a mean position £(t) and standard deviation a(t) exhibit- 
ing the same sublinear dependence on t (in the presence 
of a bias or head). A continuous time random walk 
(CTRW) transport process has successfully accounted for 
these unusual basic features. The anomalous behavior is 
an outcome of a CTRW governed by a long-tail distri- 
bution of the individual transition times or event times 
(which have to be defined in each physical context, e.g., 
trap release), 



-1-/3 



(11) 



where the Riemann-Liouville operator is defined in terms 
of the convolution [pi , 



a D]- p W{x,t) = 



1 



r(/3) dt Jo (t - t'Y-p 



o < p < 1. 

(16) 

The form of FADE ( |l5| ) is the natural generalization of 
the ADE (|sj) for power-law forms of if). The solutions of 
FADE reproduce those of the CTRW in this special case, 
and in the spatial continuum limit they form a natural 
bridge between the CTRW and the ADE [§. 



where Ap is a constant. Over the observation time range 
in which (|ll]) obtains we have for the propagating plume 
P(l,t) 



t{t) ~ t 13 , cr(t) ~t P < < 1 



(12) 



The highly dispersive nature of this P{l,t) can be dis- 
cerned by the ratio a(t)/£(t) ~ constant, instead of the 
familiar 1/yt for a Gaussian plume. The latter is pro- 
duced for > 2 (while for 1 < < 2 one has an inter- 
mediate i-dependence) . 

The CTRW accounts naturally for the cumulative ef- 
fects of a sequence of transitions which constitutes the 
particle transport. The formalism of the CTRW has been 
well documented in the literature [Berkowitz and Scher, 
1998; Metzler and Klafter, 2000]. For brevity we show 
the key equation for our purpose here in Laplace space 
(for three dimensions) 



F(s,m) = P(s,u)/P(0,u) 



s ^ (13) 



where the particle starts at the origin at t = + , P(s,m) 
is the LT of P(s,£), the probability density to find the 
particle on s at time t (plume), and F(s,u) is the LT of 
F(s,t), the probability per time for the particle to first 
arrive at s at time t. 

Explicit evaluations for F(s,u) have been determined 
for symmetric situations (e.g., plane (line) source to plane 
(line) sink) and with the use of ( |Tl] ) p0[ . The sym- 
metric cases allow a reduction of the problem to a one- 
dimensional form of F which we write as P(i; x) and 

F(u;x) = ey^{-xu p /l) (0 < < 1) (14) 

(it dimensionless) . The inverse LT of P(m;cc) has been 
expressed in terms of a class of Fox f/- functions || , but 
it is more expedient to work directly with ( |l4| ) in the 
evaluation of hit), which will be considered in the next 
section. 

Mathematically, the CTRW (in the special case using 
( |TT| ) and (0 < < 1)) is interchangeable with the time- 
fractional advection-dispersion equation (FADE) M 



dW _ x . 
-df-° Dt 







dx 



(15) 



LONG RETENTION TIMES 



We proceed to insert (|TJ) into the LT of (0) to obtain 



h(u)oc / exp(-ZM /3 /0^ = ^- /3 (l-exp(-AM /3 /0) 

(17) 

where / is the mean step distance. The expression for 
h{u) in ([It]) has been thoroughly studied in another con- 
text [Il0|; the main features are 



hir) 



t^- 1 , t <t* 

T-P- 1 , T>T* 



(18) 



(r dimensionless time) . The exponent for r > r* ensures 
that /i(r) is integrable. On a log-log plot h(r) is two 
constant slopes, 0—1, —1 — 0, with a turnover range 
between them. The center time of this range, t* , can be 
estimated as the time for the argument of the exponent 
in @ to be - 0(1) (using u - 1/r) 



H _ i 



(19) 



(The /3-factors derive from a more detailed analysis 
[Scher and Montroll, 1975, Appendix C].) For = i, 
in ( fL7j) one can do the inverse LT and obtain ?i(r) ~ 
(ttt)-2 (1 - exp(-A 2 /4/" 2 T)), clearly showing 

The change from r to t is model dependent (to be 
discussed below). We use r = v D t/le, where v is a char- 
acteristic velocity (of the velocity distribution) and £ is a 
constant. For — i, v a ~ 100 m/yr, A ~ \ km, T~ 30 m 
and e ~ \, one has (from ( |l9| ) with r* = v t* /Te), t* ~ 
10 yr, which is a reasonable time scale for the change 
from hit) ~ t^^ 1 . Hence, our scaling result ( |l8| ) agrees 
with the KFN data (|) over the measurement time range 
(> 3 decades for t< 10 yrs) with = 7 (= 1 — m). The 
i-"cutoff" for hit) in (|l8|) is algebraic, t^ 13 ^ 1 , not expo- 



nential as in (|5|)! This behavior of hit) is indicative of 
extremely long chemical retention times in catchments. 
The turnover to the t^ 13 ^ 1 dependence (in ([l^)) is a pre- 
diction of the CTRW theory and has not yet been ob- 
served. A finite value of Th depends on a change in the 
behavior of (pd|). For t ^> t* , the tail of ipit) becomes 
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steeper and for /? > 2, the intrinsic arrival time distribu- 
tion evolves to F(t; x). 

Independent tests of our transport model can be ob- 
tained with experiments using a different tracer (d) and 
boundary conditions. For injection of the tracer d at 
various catchment points {x}, the concentration in the 
stream Cd(t; l s ) is proportional to f(x— l s , t), determined 
for these boundary conditions. [This result answers the 
question raised by Stark and Stieglitz [2000] about a site- 
specific spill.] This is also a method for obtaining a value 
of (3, which determines the shape of F(x.,t). 

There are many different physical mechanisms (cf. 
Refs. above) that can give rise to the ip(t) in ( fill ) that 
has generated our results for h(t). The common features 
are the representation of the transport process as a se- 
ries of transitions [Berkowitz and Scher, 1998, 2001] and 
the encounter within this series of a sufficient number of 
transition-times that are much larger than the median 
one. These relatively few long-time events, which can be 
due to release from traps/stagnant-regions and/or pas- 
sage through low velocity segments, can have a dominant 
influence on the overall transport. The relative weight- 
ing of these events is expressed by the exponent of the 
power-law in ([□]) . 

The hydrologic environments beneath the catchment 
subregions can contain a sufficient number of these low 
velocity sections and stagnant dead-ends. The local slope 
(and interconnection 0) of the catchment side-channels 
also affects the velocity distribution. The subsurface of 
the catchment basin can be modeled as a heterogeneous 
porous medium and/or a random fracture network Q. In 
the latter the fragment-length distribution, /(s), and the 
fragment- velocity histogram, $(£), are mapped onto the 
probability rate for a transit time t (through a fragment) 
with a displacement s, 



V(s,t)oc $(£)/(*) 



(20) 



where |£| = 1/v, £ = v, t = s£ and ip(t) = X) s V'( s )i)- 
The anomalous transport observed in the fracture net- 
work particle tracking simulations is due to the power- 
law tail of <&(£) — > £~ 1 ~' 3 at large £. This is the expected 
behavior to be found in the catchment subsurface & c (0 
distribution because this hydrologic environment is the 
same as ones modeled successfully by the CTRW the- 
ory [Berkowitz and Scher, 1998, 2001]. The experiments 
outlined above can test these properties. 



CONCLUSION 

In this letter we introduce the phenomenon of anoma- 
lous transport into an analysis of the comparison of the 
power spectra of the time series of chloride in rainfall 
and in (stream) runoff. An effective distribution of travel 
times h(t) is derived (by KFN from the data) which we 



show is composed of the summation from all sites of the 
source distribution F(s, t) in the catchment. The latter, 
in turn, is composed of a sequence of transitions governed 
by a power-law (tail) distribution ([ll]) . The results agree 
naturally with the data reported from a number of catch- 
ment studies, in the sense of a power-law scaling of h(t) 
over decades of t in the observational time range (us- 
ing reasonable parameters). A consequence of our trans- 
port model is a power-law "cutoff" of h(t) predicting ex- 
tremely long chemical retention times with subsequent 
impact on the long-term effects of contaminants on these 
ecosystems. 
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